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SUMMARY 


An imolicit time-marching finite difference method is used to predict two- 
An implicit time ma & 5 number of 


An U-1.1UC w**^**£, £ / / lA 

dimensional turbulent flow at a Reynolds number of 4.4X10 

- « 4- Q -f 1 n n or OP A 


and a Mach number of 0.574 

over a shortened NACA 0012 airfoil with a trailing edge of 4.5/^ thickness and semi 
circular shape The flow is found to be unsteady but periodic in the trai ing-e g 
reeion Thus lift and drag fluctuate at small amplitudes around mean values and at 
5i®Un;t frequencies WhlU Che overnll features are in qualitative agreenent with 
tirS^eSrimental observations reported in the literature, the accuracy is limited 
^eLuse o^certain local shortcomings of the computation grid. Several recommenda- 
tlorL horL achieve improved predictions are given for future attempts. 


1. INTRODUCTION AND BACKGROUND 


Certain airfoils, helicopter blades, and modern turbine bladings often have trail- 
ing edges of a finite thickness, usually amounting to a few percent of the chor 
length^ Typically, the flow is compressible and turbulent at high Reynolds numbers, 
is calculation o^ the flow field is still an unsolved problem in the vicinity of the 
trailing edge. This is due to the limited understanding of the quite complex 

phLomeL and to the lack of powerful adequate P^'^HoLis^S^ 

In contrast to many other flow problems where an approximate prediction consistent 
iurS“Unt“ observation is achieved on the basis of the Euler equations the 
problem'^of flow over a thick trailing edge cannot be even approximately solv y 
Lancs of tnviscid theory. For subsonic flow, a stagnation point is predicted with an 
unSaliLirrreslrrpSk whereas an infinite number of solutions are found for super- 
sonic conditions. The fact that the Euler E ^ 3 ^ 

real flow behavior is probably one of the most striking (and challenging) features 

of computation of flow over a thick trailing edge. 

Therefore, a successful approach has to account for the viscous effects. Several 
calculation models for base-region flow have been developed in the ^ 

survev on several of these boundary- layer-like approaches is given in references l 
and 2. Since almost all are derived from the flow over a simple backward-facing step, 
they are based on the assumptions that 


• the separation point is known, 

• the flow reattaches to a solid wall of inclination known a priori. 
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• the wake flow is steady, and 

• the outer flow approaching the step is homogeneous. 

A major deficiency is that the interaction between the viscous-dominated 
boundary- layer and shear-layer flow, and the nearly inviscid outer flow is poorly 
accounted for. In contrast, flow over the thick trailing edge of an airfoil or 
turbine blade has the following features: 


• There are two separation points which are unknown if the base is rounded. 

• Two streams of different properties, e.g., Mach numbers, etc., merge behind 
the trailing edge; the resulting flow direction is not known a priori. 

• At high Reynolds numbers and subsonic Mach numbers, the flow is unsteady; the 
wake is of the vortex-street-type. 

• The outer flow approaching the trailing edge is not necessarily homogeneous. 

• The viscous-dominated flow near the body and in the near wake interacts 
strongly with the nearly inviscid outer flow. 


2. MOTIVATION AND DESCRIPTION OF PROBLEM 


Several problems with strong viscous-inviscid interaction and unsteady effects 
have been successfully treated by numerically solving the Reynolds-averaged Navier- 
Stokes equations or a truncated version thereof (e.g., refs. 3-6). Therefore, the 
question arises whether such methods can lead to an improved solution of the thick 
trailing-edge problem. So far, a few attempts have been made, for example: 

• Levy, Briley, and McDonald recently used an implicit technique to predict the 
laminar flow over an ellipse at a Reynolds number of 10^ and a Mach number of 
0.2. Their solution shows steady flow with two standing vortices in the near 
wake (ref. 7). 

• Fanning and Mueller solved the vorticity transport equation and Poisson's 
equation numerically and were thus able to predict the oscillating incompres- 
sible laminar wake behind the square base of a flat plate (ref. 8). 

• Waskiewicz, Shang, and Hankey recently published the solution for supersonic 
turbulent flow over the square base of a wedge-flat plate-model (ref. 9). 

They used MacCormack's fully explicit method (ref. 10). 

To the author's knowledge, the problem perhaps most often encountered in practi- 
cal application has not been tackled so far, namely subsonic compressible turbulent 
flow at high Reynolds number over an airfoil with a trailing edge of finite thickness 
and arbitrary shape. Therefore, the work described herein is dedicated to that 
specific problem. As in references 7-9, the attempt is restricted to two-dimensional 
flow. However, in contrast to references 8 and 9, it is not restricted to flat 
plates and/or square bases. 

The basic goal consists of efficiently achieving an improved prediction of the 
unsteady subsonic flow field near a thick trailing edge by using a modern numerical 
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The snecific problems connected with such a predxctxon should 
computatxon procedure. The spe attemnts and for practical engineering 

lead to recommendations for futu flow near the thick base of a 

applications hopes to 

;i:flde airitiS'lnlormatlon^or a better understanding ot the complex-eiow 

phenomena. 


3. COMPUTATION METHOD 


A large variety of conputatlon "“bo^^hased^on^the "--Ical^solutlon^of^the 

dinaterlllowrSeir application to flow fields with arbitrary boundary contours as 
they are encountered on thick round trailing edges. 

since the details of f- 

(e.g., refs. 4 and 11), only xt Revnolds-averaged Navier-Stokes equations, 

equations are the the Cartesian 

and the energy equatxon wrxtte ® ^ onto the grid line coordinates C,h. 

coordinate system ’'•y- f a mapping of the physical 

As shoun in figure 1, 'he tra nonorthogonal coordinates, onto the com- 

plane where K and n Cartesian coordinate system. It is important to 

putational plane where they form ^ ^ ooint of the airfoil at the trail- 

iote that the wake centerline, f a last P°iM of hhe a 

lag edge (see dS- 2). -a tiirsf luUoris coated twice for 

along the wake center line, ks q 11 and 12 the grid spacing in the 

reL:it:VlfertL;n uLf 

rrLt°tLrgUhe°:o:ri«:«avier-st^^^^^ 

- - Sri Fi ""°“- 

-“^pirfs'n ed^d on tha^rassure varlati™^ 

finite thickness. Unlike certarn coupling Pt°aadures, this^aPP^_^ 

Tli'ctn :i:S^t'ro:pSfmShrb™e dUflcult In unsteady flow, the application 
of the thin-layer approach Is very attractive here. 

The thin layer model Imposes two tP=trlctions 'I*® wake with 

boundary-layer-llke coordinate system must (b) the Inter- 

the n-direction (aPPta>'l™tely) normal to the flow^dlrect^ ^ 

::uanL“s“ufe7s'rf rg^'iraP m^^ "“"f uhlch S'usriirra/errrto'as a" ' 

^:r?dr:;d'wSch’'Ks1hrp™pertrth°at ^h^ 5-dlrectlon coincides with the streamwlse 
direction near the body and xn the wake. 

The Influence of turbulence Is modeled using the ar alsabralc eddy- 

vlsco^ty moLl developed by Baldwin and Loman and described in detail m 
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reference 12. As discussed in section 5, detailed measurements are not available for 
the unsteady near-wake flow. Therefore, one is forced to adopt a turbulence model 
which was not developed specifically for that type of flow; however, it is adequate 
for the region where the flow is attached to the airfoil surface. Since the Strouhal 
number of vortex shedding from thick-trailing edges is of 0 (10”^) and thus much 
lower than the dimensionless mean frequency of the large-scale turbulent eddies which 
is of 0 (10), the use of a steady or quasi-steady turbulence model is justified. 

If the same notation used in reference 4 is adopted here, the equations to be 
solved can be written as 


at 35 3n Re 3n 

where the vectors q, E, F, and S are given in detail in reference 4. They are effi- 
ciently obtained by means of the Beam-Warming delta-form approximate-factorization 
algorithm (refs. 13 and 14). It is an implicit time-marching finite difference scheme 
of second-order accuracy in space and of first- or second-order accuracy in time, 
depending on the type of differencing used to advance the solution in time. The 
results presented in section 6 are computed using the first-order time-accurate 
scheme. 


4. GRID GENERATION 


Thompson, Thames, and Mastin have developed a practical method to generate grids 
that vary smoothly and fit arbitrary boundary contours precisely (ref. 15). In its 
simplest form, the grid in the physical plane is found by solving the two Laplace 
equations 


^+^=0 

3x^ 3y^ 

3x^ 3y^ 

which, for the actual solution are transformed onto the computational plane and solved 
there by relaxation on the same grid that is used to solve the flow equations. 

This method has the advantage that the user can select the location of the grid 
points on the boundary contours (Dirichlet boundary conditions). Thus, clustering of 
the boundary nodes can easily be adjusted to the actual problem which, in the present 
case, requires a fine resolution of the trailing-edge contour. A major deficiency is 
that the mesh spacing in n-direction near the body and in the wake is far too coarse 
for viscous-flow computation. Therefore, Sorenson and Steger proposed to discard the 
grid point distribution along the lines of constant 5 and to recluster the nodes by 
a simple monotonic stretching function, which moreover allows the minimum mesh spacing 
at the surface to be specified by the user (ref. 16). According to reference 16, 
constant minimum spacing is used along the body. It may be based on a simple estima- 
tion of the boundary-layer thickness at the trailing edge by means of the well-known 
Reynolds -number criterion. In an attempt to account for the variation of the 
boundary- layer thickness along the surface, the minimum mesh spacing is continuously 
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reduced toward the front part of the airfoil. A similar modification is used in the 

wake. 


A further modification is needed because the lines of constant 5 intersect the 

=r^h: „ar:n“S 

the Laplac'^ equations. Here, however, a much simpler procedure is u • 
required that the nodes on the surface and the lines of constant 5 disnlac- 

flow" which are well distributed, are retained. The basic , 

ing the grid points near the body along the lines of constant /J* 
a straight line normal to the surface is drawn. As indicated in figu , 
poin^ if connected by another straight line with the first 

Lter part Smooth curvature is achieved by continuously cutting off the resulting 
co^Lr iwo quantities need to be specified: (a) the number of nodes which must 

belong to the strictly normal direction of the new lines of constant 5, and (b) 

here, which should be viewed as a simple engineering approach rather than a very 
sophisticated solution, avoids any iteration. 


5. CHOICE OF A TEST CASE 


Experimental data on the unsteady near wake of litSa^ur^ 

with a trailing edge of finite thickness are very rare. A search in the literatu 
reLlted in a 5ery limited collection of Strouhal numbers for flat plates and symme 
rical airfoils and is described in reference 17. 

A common feature of almost all of the experiments is that no details are avail- 
able except the Reynolds number, the Mach number, and the body shape. According o 
refLencr^rthrLequency of vortex shedding can be strongly affected by the ratio 
of boundary-layer displacement thickness to trailing-edge thickness. It is therefor 
aLon5!sS tSrmost of the experimentalists do not give any infomation on the 
boundarv layer Often, it is not even known whether the boundary layer is l^ina 
ortuSileS Therefore, we must abandon the attractive idea of confining the compu- 
tation domain to the trailing-edge region and using measured boundary layer profi e 
as upstream boundary conditions. 

Numerical computation methods give a detailed prediction of the flow quantities 
at eve^^ gSd poln?. To check these data reliably and to detect the weakness of a 
mLhirLLiLly as well as to improve these data, detailed measurements are needed 
in the Ltire region of interest. Due to the lack of such experimental data, it is 
not pLSbirto Ltablish a test case which allows reliable checks on the nt^erical 
prediction by comparing it against measurements. Therefore, ^he author prefers 
design a theoretical test case which allows several internal chec . 

For reasons of symmetry, a NACA 0012 airfoil is chosen. Its rear portion is cut 
off and replaced by a semicircular base with a radius of 2/, 

length (see fig. 4). The requirement that the trailing-edge contour fits the airt 
surfaL and its slope results in a new chord length of 87.3% of the original one 
Thus, the actual ratio of trailing-edge thickness to chord length amounts to 4.5%. 
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From the remarks just given, it follows that the entire flow field around the 
airfoil must be calculated. The resolution of the trailing-edge region requires a 
lot of grid points. The total number of grid points available is limited for obvious 
reasons . Since the resolution of shocks would require a fine grid in these regions 
and thus would reduce the number of nodes available at the base, the flow should 
preferably be subsonic everywhere. However, effects of compressibility mtist not be 
negligible. This is achieved at the free-stream Mach number of 0.574 used in the 
computation. The free-stream Reynolds number amounts to 5x10® and the boundary layer 
is assumed to be turbulent over its entire length. Two different angles of attack 
are selected: 0° and 2®. This choice allows several checks based on physical con- 

siderations, as will be seen in the discussion of the results. 

Since the governing equations are used in nondimensional fom, it should be 
noted that while the reference quantities are arbitrary, the Reynolds number is 
defined in terms of these reference values. For convenience, the chord of the origi- 
nal NACA 0012 airfoil is chosen as a reference length. The Reynolds number of 5x10® 
refers to that length. Its value based on the new chord is thus smaller by 12.7%, 
that is, it amounts to 4.36x10®. 


6 . RESULTS AND DISCUSSION 


A grid with a total number of 77 x 34 nodal points is used. The grid section 

near the body shown in figure 5 and the trailing-edge detail shown in figure 6 illus- 

trate the modification of the grid generation procedure described in section 4. 

Within the 12 nodes next to the body surface or wake center line, the lines of con- 
stant 5 are strictly normal to these contours. As is evident from figure 6 , this is 
not true for the semicircular base. There, the direction is based on a rough esti- 
mation of the mean-flow direction expected. The minimum mesh spacing in the 
n-direction increases from 0.125x10“^ at the leading edge up to 0 . 220 x 10 “® at the 
trailing edge. A check of the predicted flow field just ahead of the base indicates 
that 16 grid points are inside the boundary layer. 

Figure 6 also shows that the lines of constant n do not curve smoothly as they 

intersect the two lines of constant 5 emanating from the last nodal point at the 

trailing edge. This deficiency is likely to cause locally a certain loss of accuracy, 
a fact which was underestimated when the grid was generated. Indeed, the most 
crucial detail in creating an adequate network is the intersection of the wake center 
line with the base contour. To achieve smooth curvature, one would have to accumulate 
many grid points in that region and probably to introduce a cusp of almost vanishing 
size. Another way to avoid the deficiency is by using a so-called 0-grid where the 
lines of constant n form closed loops around the airfoil. However, it must yet be 
examined in that case whether the thin-layer approximation is still valid in the 
trailing-edge region. 

The first prediction was dedicated to symmetrical flow, that is, a = 0°. The com- 
putation was started using free-stream values at all grid points. Thus any artificial 
asymmetry was avoided. The length of the time step was doubled after 100 and 
200 steps and then kept constant. Already, after 400 time steps corresponding to 
approximately 3 - 1/2 actual chord lengths traveled, a symmetrical oscillation of the 
lift coefficient at a distinct frequency was observed. Since the amplitude and the 
frequency were still increasing, the computation was continued. After they became 
constant, several hundred additional time steps were made in order to assure that this 
was not an intermediate stage but the final solution. This is illustrated in figure 7 
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where the lift coefficient CL and the pressure drag coefficient CD are sho^ 
versus the number of actual chord lengths traveled. The analysis of the flow field 
shows that the oscillations are due to the periodic motion upwards and downwards of 
the separated flow at the base. This causes a slight alteration of the pressure dis- 
tribution up to a certain distance ahead of the trailing edge. The appearance of 
such fluctuations is in agreement with Summers and Page, who obse^ed it in their 
experiments with circular-arc airfoil sections with a blunt trailing edge (ref. ). 
Unfortunately, no pressure measurements were made which would allow a check of the 
magnitude of the oscillations. 


Since the airfoil is symmetrical and since the angle of attack is 0°, the lift 
fluctuation must be symmetrical with a mean value equal to zero. Using the 
definitions 


CL 


= i f CL(t)dt 
•^T 


CL 


max 


1 ^ 

= E CL 

n 1 


max 


CL . ^ 

min 


1 ^ 

= ^ L CL 
n " 1 


min 


where T is the length of a period and CLi„ax and Cl^m denote positive and negative 
peak values , we obtain 


CL = -0.893x10 


.-5 


CL - CL = 0.910x10 
max 

CL - ^ . = 0.909x10 

min 


.-2 


,-2 


The Strouhal number based on the trailing-edge thickness, d. 


Str = 


f X d 
U^ 


where f is the dimensional frequency, is found to be 0.070. is based the 

maximum thickness of the body, which is 13.7% due to shortening of chord, its 

value amounts to 0.24. Unlike the well-known case of a circular cylinder 
flow, the Strouhal number of slender bodies with thick trailing 

only on the Reynolds number (ref. 17). There is strong evidence that the Strouhal 
number is affected by the ratio of boundary- layer thickness to base height. is 

also likely to be affected by certain shape parameters, for example, maximum thickne 
and tail angle. As mentioned previously, boundary-layer transition is assumed at 
lading edge and thus it is fairly thick at the trailing edge. Since a reliable well- 
documented empirical correlation, such as the Strouhal number-Reynolds number rela- 
tionship of the circular cylinder, is not available for airfoils, the correctness 
the predicted Strouhal number cannot be assessed. 

According to figure 7, the pressure drag fluctuates at twice the frequency of the 
lift. This feature is due to the symmetry of the body and the choice of a 
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can also be observed on circular cylinders (see ref. 19). It expresses that the flow 
phenomena occurring in the upper part of the near wake later takes place in the lower 
half of a CL-period. Essentially, these effects are the washing- in and washing-out 
process of the separated flow and are connected with the alternating vortex shedding 
from upper and lower detachment regions. A detailed analysis of the predicted near- 
wake flow Indicates that the time-dependent pressure distribution during the washing- 
out motion differs from that during the washing- in motion. Since a substantial part 
of the pressure drag is due to the thick trailing edge, the shape of the drag fluc- 
tuation must be asymmetrical. However, two consecutive drag cycles corresponding to 
two vortices shed from different but symmetrical separation points must be identical. 
In the author's opinion, there are no physical reasons which contradict this explana- 
tion. Moreover, there is another case of self-driven unsteady flow separation with 
analogous features, namely transitory stall in diffusers. According to Kline (private 
communication) the time dependent pressure distribution during the washing-out motion 
of the separation bubble differs from the one during the washing- in process. 

Consider the same airfoil at the same Mach and Reynolds number but at a certain 
angle of attack which, however, is assumed to be small enough to prevent the separa- 
tion from moving upstream. Beside the fact that the airfoil must have a certain mean 
lift, at least the following alteration in the flow behavior can be predicted without 
any calculations: (a) the drag has to fluctuate at the same frequency as the lift, 

and (b) the Strouhal number changes only slightly. Therefore, a computation was made 
for a = 2°. All other flow parameters and program parameters (e.g., the grid, the 
time step-length, etc.) were kept unchanged. Furthermore, the computation was started 
again using free-stream values in all nodal points. 

Again, lift and drag fluctuations at a distinct frequency were observed after a 
few hundred time steps. After constant frequency and amplitudes were achieved, the 
computation of several hundred additional time steps did not change the prediction. 
This, as well as the fact that the length of drag and lift cycles are identical, is 
illustrated in figure 8. The Strouhal number based on the trailing-edge thickness 
amounts to 0.069 wl^reas the one based on the maximum thickness is 0.237. The mean 
lift coefficient CL is found to be 0.205. The values 

^ - CL = 0.7929x10"^ 

max 

CL - CL . = 0.7931x10"^ 

min 

suggest that the lift still fluctuates nearly S3mimetrically around its mean value. 

Several remarks have to be made with respect to accuracy. In order to give 
helpful suggestions for future attempts, this is done quite openly. Since the total 
number of nodes was limited, the grid points used to resolve the base contour had to 
be saved along the remaining part of the surface. This resulted in streamwise mesh 
spacings which are somewhat coarse in the mid-chord region. The solution for the last 
node at the trailing edge is computed twice independently. For example, a check on 
the numerical values indicates that the static pressures always differ somewhat; in 
the worst case, the difference is about 5%. This inaccuracy is mainly due to the 
corner-like curvature of the line of constant n near that grid point. 

Another unsatisfactory detail is illustrated in figures 9a, b, c, and d where 
four large-scale plots of the instantaneous velocity vectors in the near wake are 
shown for a = 0°. They are taken after approximately 1/12, 4/12, 7/12, and 10/12 of 
a CL-period and imply that the velocity component in the x-direction is not computed 
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correctly along the center line. This is likely at least partially due to the (too) 
large streamwise spacing in the reverse-flow region. The fact that the shortcoming 
is the worst at maximum and minimum lift (i.e., when the flow field is the most 
asymmetrical) and that it disappears at zero lift (fig. 10) may suggest that the 
boundary condition applied there is questionable. An extrapolation procedure from 
above and below followed by averaging was used in the computation. Based on these 
considerations, the accuracy of all the numerical values given above must be consid- 
ered limited. In the author's opinion, the shortcoming of the grid near the intersec- 
tion of base contour and wake center line has more negative impact than the boundary 
condition. 

In addition to the points given previously, some further recommendations for 
future attempts can be derived from the previous discussion. A finer resolution of 
the reverse-flow region is necessary. As a rule of thumb, its length may be assumed 
to be roughly equal to the base height. The stretching in the n-direction should 
preferably be reduced there (see fig. 6). The limited upstream influence of the 
unsteady flow phenomena at the trailing edge suggests that the number of grid points 
available in the near wake can be increased by confining the computation domain to the 
rear part of the airfoil. The present results imply that the inflow plane could be 
located in the mid-chord region with no adverse effects. Beside the application of 
an improved procedure for the nodes on the near-wake center line, it may also be 
useful to scrutinize the terms neglected in the thin-layer model. As mentioned by 
several authors (e.g. , in ref. 20), there is evidence that the normal stress terms in 
the streamwise direction become significant near separation. 


7. SUMMARY AND CONCLUSIONS 


An implicit time-marching finite difference scheme was used to solve the thin- 
layer equations for two-dimensional turbulent flow at a Reynolds number of 4.4x10^ 
and a Mach number of 0.574 over a modified NACA 0012 airfoil with a 4.5% thick round 
trailing edge. The computation grid was generated using an existing automatic pro- 
cedure. A simple modification was added to achieve strictly normal direction of the 
radial grid lines near the body. Solutions were computed for a = 0° and 2 °. Due to 
unsteady but periodic flow in the trailing-edge region, lift and drag fluctuate at 
distinct frequencies around mean values. The amplitudes are small compared to the 
mean values. The overall features of the numerical solution are in agreement with 
the few experimental observations known from literature. The accuracy cannot be 
assessed precisely due to the lack of detailed measurements. The accuracy is judged 
to be limited because the streamwise grid lines do not curve smoothly near the inter- 
section point of base contour and wake center line and because the number of computa- 
tion stations in the streamwise direction is small in the near wake. The recommenda- 
tions given on how to avoid these shortcomings and on how to improve some further 
details should be useful in achieving more accurate predictions in future attempts. 
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Figure 6.- Grid detail in the base region. 
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(a) After 4/50 of a CL-perlod. 



(b) After 16/50 of a CL-period. 

Figure 9.- Instantaneous velocity vectors in the base region. 
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